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Abstract. We present Path Integral Monte Carlo C code for calculation of quantum 
mechanical transition amplitudes for ID models. The SPEEDUP C code is based on the 
use of higher-order short-time effective actions and implemented to the maximal order 
p=18 in the time of propagation (Monte Carlo time step), which substantially improves 
the convergence of discretized amplitudes to their exact continuum values. Symbolic 
derivation of higher-order effective actions is implemented in SPEEDUP Mathematica 
codes, using the recursive Schrodinger equation approach. In addition to the general 
ID quantum theory, developed Mathematica codes are capable of calculating effective 
actions for specific models, for general 2D and 3D potentials, as well as for a general 
many-body theory in arbitrary number of spatial dimensions. 
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1 Introduction 

Exact solution of a given many-body model in quantum mechanics is usually expressed 
in terms of eigenvalues and eigenfunction of its Hamiltonian 

^=E^+^(qi Am), (1.1) 

but it can be also expressed through analytic solution for general transition amplitude 
A(a,b;T) = (b|e^'^^''^|a) from the initial state |a) to the final state |b) during the time of 
propagation T. Calculation of transition amplitudes is more suitable if one uses path inte- 
gral formalism fP-'S], but in principle, if eigenproblem of the Hamiltonian can be solved. 
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one should be able to calculate general transition amplitudes, and vice versa. However, 
mathematical difficulties may prevent this, and even more importantly, exact solutions 
can be found only in a very limited number of cases. Therefore, use of various analytic 
approximation techniques or numerical treatment is necessary for detailed understand- 
ing of the behavior of almost all models of interest. 

In numerical approaches it could be demanding and involved to translate numerical 
knowledge of transition amplitudes to (or from) eigenstates, but practically can be always 
achieved. It has been implemented in various setups, e.g. through extraction of the 
energy spectra from the partition function ||2}^, and using the diagonalization of space- 
discretized matrix of the evolution operator, i.e. matrix of transition amplitudes ll6]- [T0ll . 
All these applications use the imaginary-time formalism |ITn[T2| . typical for numerical 
simulations of such systems. 

Recently introduced effective action approach IfTsHlTll provides an ideal framework 
for exact numerical calculation of quantum mechanical amplitudes. It gives systematic 
short-time expansion of amplitudes for a general potential, thus allowing accurate cal- 
culation of short-time properties of quantum systems directly, as has been demonstrated 
in Refs. IISHTOl. For numerical calculations that require long times of propagation to be 
considered using e.g. Monte Carlo method, effective action approach provides improved 
discretized actions leading to the speedup in the convergence of numerically calculated 
discretized quantities to their exact continuum values. This has been also demonstrated 
in Monte Carlo calculations of energy expectation values using the improved energy es- 
timators I51IT8I. 

In this paper we present SPEEDUP codes [19J which implement the effective action 
approach, and which were used for numerical simulations in Refs. ll4ll5ll8]- [10l[T3] - [T73 . The 
paper is organized as follows. In Section |2] we briefly review the recursive approach for 
analytic derivation of higher-order effective actions. SPEEDUP Mathematica codes capa- 
ble of symbolic derivation of effective actions for a general one- and many-body theory 
as well as for specific models is described in detail in Section |3j while in Section |4] we 
describe SPEEDUP Path Integral Monte Carlo C code, developed for numerical calcula- 
tion of transition amplitudes for ID models. Section |5] summarizes presented results and 
gives outlook for further development of the code. 

2 Theoretical background 

From inception of the path integral formalism, expansion of short-time amplitudes in 
the time of propagation was used for the definition of path integrals through the time- 
discretization procedure [21131 . This is also straightforwardly implemented in the Path 
Integral Monte Carlo approaches |20 1, where one usually relies on the naive discretization 
of the action. Several improved discretized actions, mainly based on the Trotter formula 
and its generalizations, were developed and used in the past ||2T]423II . A recent analysis of 
this method can be f oimd in Jang et al Ii24| . Several related investigations dealing with the 
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speed of convergence have focused on improvements in short-time propagation [25','26] 
or the action \27\. More recently, spHt-operator method has also been developed Ii28n32j , 
later extended to include higher-order terms [33 -.36I , and systematically improved using 
the multi-product expansion Il37l - l39| . 

The effective action approach is based on the ideal discretization concept lfT6ll . It was 
introduced first for single-particle ID models l|T3l - [T5l and later extended to general many- 
body systems in arbitrary number of spatial dimensions |5, 17J. This approach allows 
systematic derivation of higher-order terms to a chosen order p in the short time of prop- 
agation. 

Recursive method for deriving discretized effective actions, established in Ref. flT], is 
based on solving the underlying Schrodinger equation for the amplitude. It has proven 
to be the most efficient tool for treatment of higher-order expansion. In this section we 
give brief overview of the recursive method, which will be implemented in Mathematica 
in the next section. We start with the case of single particle in ID, used in the SPEEDUP 
C code. Throughout the paper we will use natural system of units, in which h and all 
masses are set to unity. 



2.1 One particle in one dimension 

In the effective action approach, transition amplitudes are expressed in terms of the ideal 
discretized action S* in the form 

A{a,b;T) = -^e-'*(^''>^l (2.1) 

which can be also seen as a definition of the ideal action IIT6H . Therefore, by definition, 
the above expression is correct not only for short times of propagation, but for arbitrary 
large times T. We also introduce the ideal effective potential W, 



S*{a,b;T) = T 



1 fh-a^ ^ 



T 



-W 



(2.2) 



reminiscent of the naive discretized action, with the arguments of the effective potential 
(fl, b, T) usually written as W ^^,-^;T^ , to emphasize that we will be using mid-point 
prescription. 

However, ideal effective action and effective potential can be calculated analytically 
only for exactly solvable models, while in all other cases we have to use some approxi- 
mative method. We use expansion in the time of propagation, assuming that the time T 
is small. If this is not the case, we can always divide the propagation into N time steps, 
so that e = T/N is small. Long-time amplitude is than obtained by integrating over all 
short-time ones, 

A{a,b;T) = J dqi-'-dq^-i A{a,qi;£)A{qi,q2;£)---A{qN_i,b;£), (2.3) 
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paving the way towards Path Integral Monte Carlo calculation, which is actually imple- 
mented in the SPEEDUP C code. 

If we consider general amplitude A{q,q')£), introduce the mid-point coordinate x = 
{q+q') 1 2 and deviation x = {q' — q)/ 2, and express A using the effective potential, 

A{q,q';e) = (2.4) 

the time-dependent Schrodinger equation for the amplitude leads to the following equa- 
tion for W 

W+xdW+£d,W-led^W-lEd^W+le^{dWf + l£^{mf = l{V+ + V^), (2.5) 
8 8 8 8 2 

where y±=y(x±f), i.e. V^ = V{q), V^=V{q'). The short-time expansion assumes that we 
expand W to power series in e to a given order, and calculate the appropriate coefficients 
using Eq. (I2.5D . We could further expect that this results in coefficients depending on the 
potential V{x) and its higher derivatives. However, this scheme is not complete, since 
the effective potential depends not only on the mid-point x, but also on the deviation 
X, and the obtained equations for the coefficients cannot be solved in a closed form. In 
order to resolve this in a systematic way, we make use of the fact that, for short time of 
propagation, deviation x is on the average given by the diffusion relation x^ cx e, allowing 
double expansion of W in the form 

oo m 

W{x,x;e)=J2E''nA^)^""'^"'- (2-6) 

Restricting the above sum over m to p — 1 leads to level p effective potential Wp{x,x}e) 
which gives expansion of the effective action S* to order eP, and hence the level desig- 
nation p for both the effective action and the corresponding potential Wp. Thus, if the 
diffusion relation is applicable (which is always the case in Monte Carlo calculations), in- 
stead of the general double expansion in x and e, we are able to obtain simpler, systematic 
expansion in e only. 

As shown previously IIT3] - [T5l , when used in Path Integral Monte Carlo simulations for 
calculation of long time amplitudes according to Eq. (|2.3|l , use of level p effective action 
leads to the convergence of discretized amplitudes proportional to e^, i.e. as 1 / W, where 
N is the number of time steps used in the discretization. 

If we insert the above level p expansion of the effective potential to Eq. (|2.5|) , we obtain 
the recursion relation derived in Ref. UlTll , 

m-2 

8{m + k + l)c„,k = {2k + 2){2k+l)c,,,k+l+C_,^k-LL''ir<n-l-2,k-r 

1=0 r 

m-2 

- ^ J^2r{2k-2r+2) (2.7) 

1=1 r 
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Figure 1: Order in which the coefficients c„, j; are calculated: diagonal ones from Eq. (|2.8|1 . off-diagonal from 
recursion (|2.7p . 



where the sum over r goes from max{0, k—m + l+2} to mm{k, /}. This recursion can 
be used to calculate all coefficients Cm,k to a given level p, starting from the known initial 
condition, co,o = ^- The diagonal coefficients can be calculated immediately, 

^'"'"=(2m + l)!' ^^-^^ 

and for a given value of m = 0, . . . p — 1, the coefficients c,„ jt follow recursively from evalu- 
ating i2.7^ for = m — 1, . . . , 1,0, as illustrated in Fig. [TJ 



2.2 Extension to many-body systems 

The above outlined approach can be straightforwardly applied to many-body systems. 
Again the amplitude is expressed through the effective action and the corresponding 
effective potential, which now depends on mid-point positions and deviations of all par- 
ticles. For simplicity, these vectors are usually combined into D x M dimensional vectors 
X and X, where D is spatial dimensionality, and M is the number of particles. In this 
notation, 

A(qq'-£) = - g-7^^-fW(w) (29) 

where initial and final position q= (qi,...,qM) and q' = (q'i,---/qA4) ^re analogously de- 
fined D X M dimensional vectors. Here we will not consider quantum statistics of parti- 
cles. The required symmetrization or antisymmetrization must be applied after transition 
amplitudes are calculated using the effective potential. 

Many-body transition amplitudes satisfy D x M-dimensional generalization of the 
time-dependent Schrodinger equation, which leads to the equation for the effective po- 
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tential similar to Eq. I|2.5|) . with vectors replacing previously scalar quantities, 

1 1 T 1 1 

w+x-8w+£a,w--£a%--£a^w+-£2(aw)2+-£2(aw)2=-(y++y_). (2.10) 
000 2 

The effective potential for short-time amplitudes again can be written in the form of the 
double expansion in e and x. However, it turns out to be advantageous to use the expan- 
sion 

00 m 

W{x,x;e)= J2 J2^"'-'W,„,k{x,x) , (2.11) 

m=Qk=0 

and work with fully contracted quantities W,„ )t 

W„,,, ) = x,^ x,^ ■ ■ ■ x,J^^^{x) , (2.12) 

rather than with the respective coefficients c^'^ ''^^ In this way we avoid the computa- 
tionally expensive symmetrization over all indices ii,...,i2k- After inserting the above 
expansion into the equation for the effective potential, we obtain the recursion relation 
which represents a generalization of previously derived Eq. (12.711 for ID case, and has the 
form 

m-2 

8{m+k+l)W^,k = d^W^-u + rW„r,k+i-J2J2{dWi,r)-{dW^-i-2,k-r) 

- E E(^^'.'-)-(^^— '-U-r+i)- (2-13) 

1=1 r 

The sum over r runs from max{0, k—m + l+2} to mm{k, I}, while diagonal quantities 
Wm,m can be calculated directly, 

w.,.=^^^(x.a)^'«y. (2.14) 

The above recursion disentangles, in complete analogy with the previously outlined case 
of one particle in ID, and is solved in the order shown in Fig.[l] 



3 SPEEDUP Mathematica codes for deriving the higher-order 
effective actions 

The effective action approach can be used for numerically exact calculation of short-time 
amplitudes if the effective potential Wp can be analytically derived to sufficiently high 
values of p such that the associated error is smaller than the required numerical pre- 
cision. The error eP for the effective action, obtained when level p effective potential 
is used, translates into eP^C)M/2 £qj. ^ general many-body short-time amplitude. How- 
ever, when amplitudes are calculated using the Path Integral Monte Carlo SPEEDUP C 
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code [m, which will be presented in the next section, the errors of numerically calculated 
amplitudes are always proportional to eP ~ 1/N^, where N is number of time-steps in the 
discretization of the propagation time T. 

Therefore, accessibility of higher-order effective actions is central to the application of 
this approach if it is used for direct calculation of short-time amplitudes ||8HT0ll, as well 
as in the case when PIMC code is used ll4ll51 [T8ll . However, increase in the level p leads 
to the increase in complexity of analytic expressions for the effective potential. On one 
hand, this limits the maximal accessible level p by the amount of memory required for 
symbolic derivation of the effective potential. On the other hand, practical use of large 
expressions for Wp may slow down numerical calculations, and one can opt to use lower 
than the maximal available level p when optimizing total CPU time required for numer- 
ical simulation. The suggested approach is to study time-complexity of the algorithm in 
practical applications, and to choose optimal level p by minimizing the execution time 
required to achieve fixed numerical precision. 

We have implemented efficient symbolic derivation of higher-order effective actions 
in Mathematica using the recursive approach. All source files described in this section 
are located in the Mathematica directory of the SPEEDUP code distribution. 

3.1 General ID Mathematica code 

SPEEDUP code [19J for symbolic derivation of the effective potential to specified level 
p is implemented in Mathematica [40], and is available in the EffectiveAction-lD.nb 
notebook. It implements the algorithm depicted in Fig. [T] and calculates the coefficients 
c„,j( for ra = 0,...,p — 1 and k = m,...,0, starting from the initial condition co,o = ^- For a 
given value of tn, the diagonal coefficient Cttj^m 

is first calculated from Eq. (I2.8D , and then 
all off-diagonal coefficients are calculated from the recursion (|2.7|) . 

In this code the potential V{x) is not specified, and the effective potential is derived 
for a general one-particle ID theory. The resulting coefficients j. and the effective po- 
tential are expressed in terms of the potential V and its higher derivatives. Level p effec- 
tive potential, constructed as 

Wp{x,x;e) = Y^c,,^{xy^-H'-\ (3.1) 

m=Ofc=0 

contains derivatives of V to order 2p — 2. 

The only input parameter of this Mathematica code is the level p to which the effec- 
tive potential should be calculated. As the code runs, it prints used amount of memory 
(in MB) and CPU time. This information can be used to estimate the required comput- 
ing resources for higher values of p. The calculated coefficients can be exported to a 
file, and later imported for further numerical calculations. As an illustration, the file 
Ef f ectiveActioii-lD-export-p5 .m contains exported definition of all the coefficients 
c„, /f calculated at level p = 5, while the notebook Ef f ectiveAction-lD-matching-p5 . nb 
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contains matching output from the interactive session used to produce the above p = 5 
result. 

The execution of this code on a typical 2 GHz CPU for level p = 10 requires 10-15 MB 
of RAM and several seconds of CPU time. We have successfully run this code for levels as 
high as p = 35 |fT9i . SPEEDUP C code implements effective actions to the maximal level 
p = 18, with the size of the corresponding C function around 2 MB. If needed, higher 
levels p can be easily implemented in C and added to the existing SPEEDUP code. 

3.2 General 2D and 3D Mathematica code 

Although we have developed Mathematica code capable of deriving effective actions 
for a general many-body theory in arbitrary number of spatial dimensions, in practical 
applications in 2D and 3D it can be very advantageous to use simpler codes, able to 
produce results to higher levels p than the general code |[9l[T0|. 

This is done in files Eff ectiveAction-2D.nb and EffectiveAction-3D.nb, where 
the recursive approach is implemented directly in 2D and 3D. Execution of these codes 
requires more memory: for p = 10 effective action one needs 60 MB in 2D case, while 
in 3D case the needed amount of memory increases to 860 MB. On the other hand, the 
execution time is several minutes for 2D code and around 30 minutes for 3D code. 

The distribution of the SPEEDUP code contains exported p = 5 definitions of con- 
tractions W„, )t for both 2D and 3D general potential, as well as matching outputs from 
interactive sessions used to generate these results. 

3.3 Model-specific Mathematica codes 

When general expressions for the effective actions, obtained using the above described 
SPEEDUP Mathematica codes, are used in numerical simulations, one has to specify the 
potential V and its higher derivatives to order 2p — 2 in order to be able to calculate tran- 
sition amplitudes. Such approach is justified for systems where the complexity of higher 
derivatives increases. However, for systems where this is not the case, or where only a 
limited number of derivatives is non-trivial (e.g. polynomial interactions), it might be 
substantially beneficial to specify the potential at the beginning of the Mathematica code 
and calculate the derivatives explicitly when iterating the recursion. 

Using this approach, one is able to obtain coefficients c„, and the effective poten- 
tial W directly as functions of the mid-point x. This is implemented in the notebooks 
Ef f ectiveAction-lD-AHO .nb and Ef f ectiveAction-2D-AH0 .nb for the case of anhar- 
monic oscillators in ID and 2D, 



^2D-AHo{x) 



^lD-AHo(^) 




(3.2) 



(3.3) 
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These codes can be easily executed within few seconds and with the minimal amounts of 
memory even for p = 20. For ID anharmonic oscillator we have successfully calculated 
effective actions to excessively large value p = 144, and in 2D to p = 67 I19J , to illustrate 
the advantage of this model-specific method. 

Similar approach can be also used in another extreme case, when the complexity of 
higher derivatives of the potential V increases very fast, so that entering the correspond- 
ing expressions to the code becomes impractical. Even in this situation expressions for 
effective actions can be usually simplified using some appropriate model-specific ansatz. 
The form of such ansatz can be deduced from the form of model-specific effective po- 
tentials, and then used to simplify their derivation. Such use-case is illustrated in the 
SPEEDUP Mathematica code for the modified Poscl-Teller potential. 



A 

{coshoix) 



Vid-mpt{x) = - 2 - (3.4) 



For this potential, the coefficients Cjn,k of the effective potential can be expressed in the 
form 

, , ^ J (tanhax)^' 
c...(x) = Erf„a; (^^^t^^^)2.-2/+2 - (3.5) 

and newly introduced constant coefficients dfn,k,i can be calculated using the model-speci- 
fic recursion in Ef f ectiveAction-lD-MPT.nb. The form of the ansatz (|3.5D is deduced 
from the results of executing general ID Mathematica code, with the model-specific po- 
tential (13.41) defined before the recursion calculation of the coefficients is performed. Us- 
ing this approach, we were able to obtain maximal level /3 = 41 effective action |1T9| . 



3.4 General many-body Mathematica code 

SPEEDUP Mathematica code for calculation of effective action for a general many-body 
theory is implemented using the MathTensor BDl package for tensorial calculations in 
Mathematica. This general implementation required some new functions related to the 
tensor calculus to be defined in the source notebook Ef f ectiveAction-ManyBody . nb pro- 
vided with the SPEEDUP code. 

The function GenNewInd [n] generates the required number n of upper and lower in- 
dices using the MathTensor function UpLo, with the assigned names upl, lol, . . . , as well 
as lists upi and loi, each containing n strings corresponding to the names of generated 
indices. These new indices are used in the implementation of the recursion for calcula- 
tion of derivatives of W,„ jt, contractions of the effective potential, and for this reason had 
to be explicitly named and properly introduced. 

The expressions obtained by iterating the recursion contain large numbers of contrac- 
tions, and function NewDef Unique [contr] replaces all contracted indices with the newly- 
introduced dummy ones in the contraction contr, so that they do not interfere with the 
calculation of derivatives in the recursion. This is necessary since the derivatives in re- 
cursion do not distinguish contracted indices from non-contracted ones if their names 
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happen to be generated by the function GenNewInd. Note that the expression contr does 
not have to be full contraction, i.e. function NewDef Unique will successfully act on ten- 
sors of any kind if they have contracted indices, while it will leave them unchanged if no 
contractions are present. 

The function NewDerivativeVec [contr , vec, ind] implements calculation of the 
first derivative of the tensor contr (which may or may not contain contracted indices, 
but if it does, they are supposed to be uniquely defined dummy ones, which is achieved 
using the function NewDef Unique). The derivative is calculated with respect to vector vec 
with the vectorial index ind. The index ind can be either lower or upper one, and has to 
be generated previously by the function GenNewInd. 

Finally, the function NewLaplacianVec [contr , vec] implements the Laplacian of the 
tensor contr with respect to the vector vec, i.e. it performs the calculation of contractions 
of the type 

r contr. (3.6) 

ovecj dvec' 

After all described functions are defined, the execution of the code proceeds by set- 
ting the desired level of the effective action p, generating the needed number of named 
indices using the function call GenNewInd [2 p + 2] , and then by performing the recur- 
sion according to the scheme illustrated in Fig.[ll The use of MathTensor function CanAll 
in the recursion ensures that the obtained expressions for W [m , k] will be simplified if 
possible. This is achieved in MathTensor by sorting and renaming all dummy indices 
using the same algorithm and trying to simplify the expression obtained in such way. By 
default, Mathematica will distinguish contracted indices in two expressions if they are 
named differently, and MathTensor works around it using the renaming scheme imple- 
mented in CanAll. 

The computing resources required for the execution of the many-body SPEEDUP 
Mathematica code depend strongly on the level of the effective action. For example, 
for level p = 5 the code can be run within few seconds with the minimal memory re- 
quirements. The notebook with the matching output of this calculation is available as 
Ef f ectiveAction-ManyBody-matching-p5 .nb, and the exported results for W[m, k] are 
available in Ef f ectiveAction-ManyBody-export-p5 .m. We were able to achieve maxi- 
mal level p = 10 HH, with the CPU time of around 2 days on a recent 2 GHz processor. 
The memory used by Mathematica was approximately 1.6 GB. 

Note that exporting the definition of the effective potential from Mathematica to a file 
will yield lower and upper indices named 111, uul, etc. In order to import previous re- 
sults and use them for further calculations with the provided Mathematica code, it is nec- 
essary to replace indices in the exported file to the proper index names used by the func- 
tion GenNewInd. This is easily done using find/ replace feature of any text editor. Prior to 
importing definition of the effective potential, it is necessary to initialize MathTensor and 
all additional functions defined in the notebook Ef f ectiveAction-ManyBody .nb, and to 
generate the needed number of named indices using the function call GenNewInd [2p+2] . 
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4 SPEEDUP C codes for Monte Carlo calculation of ID transi- 
tion amplitudes 



For short times of propagation, the effective actions derived using the above described 
Mathematica codes can be directly used. This has been extensively used in Refs. OHl, 
where SPEEDUP codes were applied for numerical studies of several lower-dimensional 
models and calculation of large number of energy eigenvalues and eigenf unctions. The 
similar approach is used in Ref . 1 10 1, where SPEEDUP code was used to study properties 
of fast-rotating Bose-Einstein condensates in anharmonic trapping potentials. The avail- 
ability of a large number of eigenstates allowed not only precise calculation of global 
properties of the condensate (such as condensation temperature and ground state occu- 
pancy), but also study of density profiles and construction of time-of-flight absorption 
graphs, with the exact quantum treatment of all available eigenf unctions. 

However, in majority of applications the time of propagation cannot be assumed to 
be small. The effective actions are found to have finite radius of convergence [8], and if 
the typical propagation times in the considered case exceed this critical value. Path Inte- 
gral Monte Carlo approach must be used in order to accurately calculate the transition 
amplitudes and the corresponding expectation values iHUm. As outlined earlier, in this 
case the time of propagation T is divided into N time steps, such that £ = T /N is suffi- 
ciently small and that the effective action approach can be used. The discretization of the 
propagation time leads to the following expression for the discretized amplitude 



and qo = a,qN = b, = {qk+i +qk) /2, Xj, = {qk+i - qk) /2. 

Level p discretized effective action is constructed from the corresponding effective 
potential Wp, calculated as power series expansion to order e^^^. Since it enters the action 
multiplied by e, this leads to discretized actions correct to order eP, i.e. with the errors of 

the order eP+^. The long-time transition amplitude A^j^\a,b;T) is a product of N short- 
time amplitudes, and its errors are expected to scale as N-eP^^^l/ , as has been shown 
in Refs. ISllTSHTSI for transition amplitudes, and in Refs. [S'.TS] for expectation values, 
calculated using the corresponding consistently improved estimators. 

4.1 Algorithm and structure of the code 

SPEEDUP C source is located in the src directory of the code distribution IIT9ll . It uses the 
standard Path Integral Monte Carlo algorithm for calculation of transition amplitudes. 




(4.1) 



where S]^ stands for the discretized level p effective action. 




(4.2) 
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The trajectories are generates by the bisection algorithm ^20], hence the number of time- 
steps N is always given as a power of two, N=2\ When the amplitude is calculated with 
2^ time steps, we can also easily calculate all discretized amplitudes in the hierarchy 2^^\ 
. . . , 2*^ at no extra cost. This requires only minor additional CPU time and memory, since 
the needed trajectories are already generated as subsets of maximal trajectories with 2^ 
time-steps. 

The trajectory is constructed starting from bisection level n = 0, where we only have 
initial and final position of the particle. At bisection level n = l the propagation is divided 
into two time-steps, and we have to generate coordinate q of the particle at the moment 
T/2, thus constructing the piecewise trajectory connecting points a at the time t = 0,q 
at t = T/2, and b at t = T. The coordinate q is generated from the Gaussian probability 
density function centered at {a + b)/2 and with the width ci = vT/2. The procedure 
continues iteratively, and each time a set of points is added to the piecewise trajectory. At 
each bisection level n the coordinates are generated from the Gaussian centered at mid- 
point of coordinates generated at level n — 1, with the width t7„ = \/T/2". To generate 
numbers // from the Gaussian centered at zero we use Box-Miiller method. 



?/ = ^-2c7-2ln^icos27r^2, (4.3) 

where numbers and ^2 are generated from the uniform distribution on the interval 
[0,1], using the SPRNG library Il42l . If the target bisection level is s, then at bisection 
level n <s we generate 2"^^ numbers using the above formula, and construct the new 
trajectory by adding to already existing points the new ones, according to 

,[(i+2o-2«-»i=,,+it^!::M±lH:!:!l. (4.4) 

where / runs from to 2"-!-!. This ensures that at bisection level s we get trajectory 
with N = 2^ time-steps, consisting of N+1 points, with boundary conditions ^j[0] =fl and 
q[N] = b. At each lower bisection level n, the trajectory consists of 2" + l points obtained 
from the maximal one (level s trajectory) as a subset of points q[i-2^^"] for z = 0,1,...,2". 

The use of trajectories generated by the bisection algorithm requires normalization 
factors from all Gaussian probability density functions with different widths to be taken 
into account. This normalization is different for each bisection level, but can be calculated 
easily during the initialization phase. 

The basic C code is organized in three source files, main.c, p.c and potential, c, 
with the accompanying header files. The file potential . c (its name can be changed, and 
specified at compile time) must contain a user-supplied function VO () , defining the po- 
tential V. For a given input value of the coordinate, VO () should initialize appropriate 
variables to the value of the potential V and its higher derivatives to the required or- 
der 2p — 2. When this file is prepared, SPEEDUP code can be compiled and used. The 
distributed source contains definition of ID-AHO potential in the file potential . c, the 
same as in the file ID-AHO . c. 
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The execution of the SPEEDUP code starts with the initialization and allocation of 
memory in the mainO function, and then the array of amplitudes and associated MC 
error estimates for each bisection level n = 0,...,s is calculated by calling the function 
mc(). After printing the output, mainO deallocates used memory and exits. Function 
mc ( ) which implements the described MC algorithm is also located in the file main . c, as 
well as the function distrO, which generates maximal (level s) trajectories. 

The function mc() contains main MC sampling loop. In each step new level s trajec- 
tory is generated by calling the function distr () . Afterwards, for each bisection level n, 
function f unc () is invoked. This function is located in the file p . c, and returns the value 
of the function e^^ , properly normalized, as described earlier. This value (and its square) 
is accumulated in the MC loop for each bisection level n and later averaged to obtain the 
estimate of the corresponding discretized amplitude and the associated MC error. 

The function f unc ( ) makes use of C implementation of earlier derived effective ac- 
tions for a general ID potential. For a given trajectory at the bisection level n, funcO 
will first initialize appropriate variables with the values of the potential and its higher 
derivatives (to the required level 2p — 2) by calling the user-supplied function VO ( ) , lo- 
cated in the file potential . c. Afterwards the effective action is calculated according to 
Eq. (I4.2D , where the effective potential is calculated by the function Wp ( ) , located in the 
file p . c. The desired level p of the effective action is selected by defining the appropriate 
pre-processor variable when the code is compiled. 

In addition to this basic mode, when SPEEDUP code uses general expression for level 
p effective action, we have also implemented model-specific mode, described earlier. If 
effective actions are derived for a specific model, then user can specify an alternative p . c 
file to be used within the directory src/models/<model>, where <model> corresponds 
to the name of the model. If this mode is selected at compile time, the compiler will 
ignore p . c from the top src directory, and use the model-specific one, defined by the user. 
The distributed source contains model definitions for ID-AHO and ID-MPT potentials 
in directories src/models/lD-AHO and src/models/lD-MPT. Note that in this mode the 
potential is specified directly in the definition of the effective potential, and therefore the 
function VO ( ) is not used (nor the potent ial . c file). 

4.2 Compiling and using SPEEDUP C code 

SPEEDUP C source can be easily compiled using the Makef ile provided in the top direc- 
tory of the distribution. The compilation has been thoroughly tested with GNU, Intel and 
IBM XLC compilers. In order to compile the code one has to specify the compiler which 
will be used in the Makefile by setting appropriately the variable COMPILER, and then to 
proceed with the standard command of the type make <target>, where <target> could 
be one of all, speedup, sprng, clean-all, clean-speedup, clean-sprng. 

The SPRNG library [42J is an external dependency, and for this reason it is located in 
the directory src/deps/sprng4 . 0. In principle, it has to be compiled only once, after the 
compiler has been set. This is achieved by executing the command make sprng. After- 
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wards the SPEEDUP code can be compiled and easily linked with the already compiled 
SPRNG library. Note that if the compiler is changed, SPRNG library has to be recompiled 
with the same compiler in order to be successfully linked with the SPEEDUP code. 

To compile the code with level p=10 effective action and user-supplied function VO () 
located in the file src/lD-AHO . c, the following command can be used: 

make speedup P=10 P0TENTIAL=1D-AH0 . c 
If not specified, PDTENTIAL=potential . c is used, while the default level of the effective 
action is P=l. To compile the code using a model-specific definition of the effective poten- 
tial, instead of the POTENTIAL variable, we have to appropriately set the MODEL variable 
on the command line. For example, to compile the supplied p . c file for ID-MPT model 
located in the directory src/models/lD-MPT using the level p = 5 effective action, the fol- 
lowing command can be used: 

make speedup P=5 M0DEL=1D-MPT 
All binaries compiled using the POTENTIAL mode are stored in the bin directory, while the 
binaries for the MODEL mode are stored in the appropriate bin/models/<model> directory. 
This information is provided by the make command after each successful compilation is 
done. 

The compilation is documented in more details in the supplied README . txt files. The 
distribution of the SPEEDUP code also contains examples of compilation with the GNU, 
Intel and IBM XLC compilers, as well as matching outputs and results of the execution 
for each tested compiler, each model, and for a range of levels of the effective action p. 
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Figure 2: Convergence of SPEEDUP Monte-Carlo results for the transition amplitude ( — 0.5,0.5;!) of ID- 
MPT potential as a function of the number of time steps N, calculated with level p = l,2,10 effective actions, 
with the parameters of the potential \ = ol = 1. The full lines give the fitted functions (|4.5|) . where the constant 
term Ap corresponds to the continuum-theory amplitude A(— 0.5,0.5;1). The number of Monte-Carlo samples 

was Nmc = 10^- 
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Once compiled, the SPEEDUP code can be used to calculate long-time amplitudes of 
a system in the specified potential V. If executed without any command-line arguments, 
the binary will print help message, with details of the usage. The obligatory arguments 
are time of propagation T, initial and final position a and b, maximal bisection level s, 
number of MC samples Nmc and seed for initialization of the SPRNG random number 
generator. All further arguments are converted to numbers of the double type and made 
available in the array par to the function V0(), or to the model-specific functions in the 
file src/iiiodels/<model>/p . c. The output of the execution contains calculated value of 
the amplitude for each bisection level n = 0,...,s and the corresponding MC estimate of 
its error (standard deviation). At bisection level n = 0, where no integrals are actually 
calculated and the discretized N= 1 amplitude is simply given by an analytic expression, 
zero is printed as the error estimate. 

Fig. |2] illustrates the typical results obtained from the SPEEDUP code on the example 
of ID-MPT theory. In this figure we can see the convergence of numerically calculated 
amplitudes with the number of time-steps N to the exact continuum value, obtained in 
the limit N — )■ oo. Such convergence is obtained for each level p of the effective action 
used. However, the convergence is much faster when higher-order effective action is 
used. Note that all results corresponding to the one value of level p on the graph are 
obtained from a single run of the SPEEDUP code with the maximal bisection level s = 10. 
The simplest way to estimate the continuum value of the amplitude is to fit numerical 
results from single run of the code to the appropriate level p fitting function IfTsHTSlI , 
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Figure 3: (left) The anharmonic potential ID-AHO, its energy eigenvalues (horizontal lines) and eigenfunctions, 
obtained by direct diagonalization of the space-discretized matrix of the evolution operator with level p = 21 
effective action and parameters A = l, §- = 48. The discretization cutoff was L = 8, spacing A = 9.76- 10^*, and 
time of propagation t=0.02. (right) Results for the double-well potential, A = -10, g=12, L = 10, A = 1.22-10"3^ 
t = 0.1. On both graphs, left y-axis corresponds to V{x) and energy eigenvalues, while scale on the right y-axis 
corresponds to values of eigenfunctions, each vertically shifted to level with the appropriate eigenvalue. 
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The constant term obtained by fitting corresponds to the best estimate of the exact ampli- 
tude which can be found from the available numerical results. 

As mentioned earlier, the effective action approach can be used for accurate calcula- 
tion of a large number of energy eigenstates and eigenvalues by diagonalization of the 
space-discretized matrix of transition amplitudes ||6tiT0|. Fig.|3]illustrates this for the case 
of an anharmonic and double-well potential. The graph on the left gives several eigenval- 
ues and eigenstates for ID-AHO potential with A = l and quartic anharmonicity g = 48, 
while the graph on the right gives low-lying spectrum and eigenfunctions of the double- 
well potential, obtained for A = —10, with the moderate anharmonicity g = 12. More 
details on this approach, including study of all errors associated with the discretization 
process, can be found in Refs. I8,i9j. 



5 Conclusions 

In this paper we have presented SPEEDUP Mathematica and C codes, which implement 
the effective action approach for calculation of quantum mechanical transition ampli- 
tudes. The developed Mathematica codes provide an efficient tool for symbolic deriva- 
tion of effective actions to high orders for specific models, for a general ID, 2D and 3D 
single-particle theory, as well as for a general many-body systems in arbitrary number of 
spatial dimensions. The recursive implementation of the code allows symbolic calcula- 
tion of extremely high levels of effective actions, required for high-precision calculation 
of transition amplitudes. 

For calculation of long-time amplitudes we have developed SPEEDUP C Path Integral 
Monte Carlo code. The C implementation of a general ID effective action to maximal 
level p = 18 and model-specific effective actions provide fast 1/NP convergence to the 
exact continuum amplitudes. 

Further development of the SPEEDUP C codes will include parallelization using MPI, 
OPENMP and hybrid programming model, C implementation of the effective potential 
to higher levels p, as well as providing model-specific effective actions for relevant po- 
tentials, including many-body systems. 
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